Dynamics and stability of vortex- ant i vortex fronts in type II superconductors 
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The dynamics of vortices in type II superconductors exhibit a variety of patterns whose origin is 
poorly understood. This is partly due to the nonlinearity of the vortex mobility which gives rise 
to singular behavior in the vortex densities. Such singular behavior complicates the application of 
standard linear stability analysis. In this paper, as a first step towards dealing with these dynamical 
phenomena, we analyze the dynamical stability of a front between vortices and antivortices. In 
particular we focus on the question of whether an instability of the vortex front can occur in the 
absence of a coupling to the temperature. Borrowing ideas developed for singular bacterial growth 
fronts, we perform an explicit linear stability analysis which shows that, for sufficiently large front 
velocities and in the absence of coupling to the temperature, such vortex fronts are stable even 
in the presence of in-plane anisotropy. This result differs from previous conclusions drawn on the 
basis of approximate calculations for stationary fronts. As our method extends to more complicated 
models, which could include coupling to the temperature or to other fields, it provides the basis for 
a more systematic stability analysis of nonlinear vortex front dynamics. 

PACS numbers: 74.25.Qt,05.45.-a 



I. INTRODUCTION 



Motivation 



The properties of type II superconductors have been 
studied extensively in past decades. The analysis of pat- 
terns in the magnetic flux distribution has generally fo- 
cused on equilibrium vortex phases. The interplay of 
pinning and fluctuation effects, especially in the high- Tc 
superconductors, gives rise to a rich variety of phases 
whose main features are by now rather well understood 
Q, Q. In comparison with equilibrium behavior, how- 
ever, our understanding of the dynamics of vortices, and 
the dynamical formation of vortex patterns, is still much 
less well developed. 

Recently, experiments with magneto-optical tech- 
niques on flux penetration in thin films have revealed the 
formation of a wide variety of instabilities. An example is 
the nucleation of dendrite-like patterns in Nb and MgB2 
films IE ^ IS ' These complex structures consist of alter- 
nating low and high vortex density regions and are found 
in a certain temperature window. Likewise, flux penetra- 
tion in the form of droplets separating areas of different 
densities of vortices has been observed in NbSe2 . Pat- 
terns with branchlike structures have been found also in 
high- Tc materials, like YBa2Cu30'!-x 0- In addition, 
the scaling of the fluctuations of a (stable) vortex front 
penetrating a thin sample has been studied 

Usually the occurrence of dendrite-like patterns in 
interfacial growth phenomena can be attributed to a 
diffusion-driven, long-wavelength instability of a straight 
front, similar to the MuUins-Sekerka instability found 
in crystal growth. In this paper we therefore investigate 
the stability of a straight front of vortices and antivortices 
which propagate into a type II superconductor. Further- 



more, according to the experimental data [TollTlL ll^ IT^ . 

the boundary between vortices and antivortices exhibits 
many features suggestive of a long-wavelength instability. 

The nucleation of dendrites associated with the prop- 
agation of a flux front into a virgin sample has been at- 
tributed to such an interfacial instability. This results 
from a thermo-magnetic coupling ^ ^ .14. .15.J where a 
higher temperature leads to a higher mobility, enhanced 
flux flow, and hence a larger heat generation. However, 
the cause of the instability at the boundary between 
fluxes of op pos ite sign is still being debated. Shapiro and 
co-workers |16| attribute these patterns to a coupling to 
the temperature fleld via the heat generated by the anni- 
hilation of vortices with antivortices. On the other hand. 
Fisher et al. 0, claim that an in-plane anisotropy of 
the vortex mobility is sufficient to generate an instability. 

There are several reasons to carefully reinvestigate the 
idea of an anisotropy-induced instability of propagating 
vortex-antivortex fronts. First of all, even though this 
mechanism was claimed to be relevant for the "turbu- 
lent" behavior at the boundaries of opposite flux regions, 
the critical anisotropy coefficients found on the basis of an 
approximation '17*, IS] correspond to an anisotropy too 
high to describe a realistic situation, even when a non- 
linear relation between the current and the electric field 
was considered. Eilllll. Secondly, the calculation 
was effectively done for a symmetric stationary interface, 
rather than a moving one. Thirdly, the physical picture 
that has been advanced |0| for the anisotropy-induced 
instability is that of a shear-induced Kelvin-Hclmholtz 
instability, familiar from the theory of fluid interfaces 
|2^. However, it is not clear how far the analogy with 
the Kelvin-Helmholtz instability actually extends. 

In order to try to settle the mechanism that underlines 
such phenomena, we investigate here the linear stability 
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of the interface between vortices and antivortices with- 
out any approximations in the case where the front of 
vortices propagates with a finite velocity. We perform 
an exphcit hnear stabihty analysis which shows that, in 
the presence of an in-plane anisotropy, vortex fronts with 
sufficiently large speed are stable in the absence of cou- 
pling to the temperature. We shall see that the issue of 
the stability of fronts between vortices and antivortices 
is surprisingly subtle and rich: while we confirm the find- 
ing of Fisher et al. 0, 0| that stationary fronts have 
an instability to a modulated state, our moving fronts 
are found to be stable for all anisotropics. Moreover, our 
calculations indicate that the stability of such fronts de- 
pends very sensitively on the distribution of antivortices 
in the domain into which the front propagates, so it is 
difficult to draw general conclusions. 

Besides the intrinsic motivation to understand this 
anisotropy issue, there is a second important motivation 
for this work. Our coarse-grained dynamics of the vor- 
tex densities is reminiscent of reaction-diffusion equations 
with nonlinear diffusion. This makes the coarse-grained 
vortex dynamics very different from the Gaussian diffu- 
sive dynamics of a linear diffusion equation. For exam- 
ple, the fact that vortices penetrate a sample with linear 
density profiles is an immediate consequence of this. 
More fundamentally, the dynamically relevant fronts in 
such equations with nonlinear diffusion are usually asso- 
ciated with nonanalytic (singular) behavior of the vor- 
tex densities — such singular behavior has been stud- 
ied in depth for the so-called porous medium equation 
I25I f2fil | , which has a similar nonlinear diffusion. In 
the case we will study, the front corresponds to a line on 
one side of which one of the vortex densities is nonzero, 
while on the other side it vanishes identically. In the 
regime on which we will concentrate, this vortex density 
vanishes linearly near the singular line. But for other 
cases encountered in the literature 0, even more 
complicated non-linear dynamical equations arise that 
are reminiscent of reaction-diffusion type models in other 

tsical systems. The case of bacterial growth models 
illustrates that the non-linearity of the diffusion 
process can have a dramatic effect on the front stability, 
so a careful analysis is called for. Nevertheless, in our 
case nonlinear diffusion by itself does not lead to an in- 
stability of the front, unlike in the bacterial growth case 
[29| or viscous fingering j^] . 

From a broader perspective, we see this work as a 
first step towards a systematic analysis of moving vor- 
tex fronts. The linear stability analysis which we will 
develop can equally well be applied to dynamical models 
which include coupling to the temperature or in which 
the current-voltage characteristic is nonlinear. For this 
reason, we present the analysis in some detail for the rela- 
tively simple case where the vortex velocity is linear with 
respect to the magnetic field gradient and the current. 
Even then, as we shall see, the basic uniformly translat- 
ing front solutions can still have surprisingly complicated 
behavior. We find that the density of vortices which pen- 



etrate the sample vanishes linearly for large enough front 
velocities, but with a fractional exponent for front veloc- 
ities below some threshold velocity [s^l- Since the latter 
regime appears to be physically less relevant, and since 
we do not want to overburden the paper with mathemat- 
ical technicalities, we will focus our analysis on the first 
regime. As stated before, in this regime we find that an 
anisotropy in the mobility without coupling to the tem- 
perature does not give rise to an instability of the flux 
fronts. 

Our analysis will be aimed at performing the full sta- 
bility analysis of the fronts in the coupled continuum 
equations for the vortex densities. Our procedure thus 
differs from the one of m lisj in which a sharp inter- 
face limit was used. In many physical systems it is of- 
ten advantageous to map the equations onto a moving 
boundary effective interface problem, in which the width 
of the transition zone for the fields is neglected. One can 
in principle derive the proper moving boundary approx- 
imation from the continuum equations with the aid of 
singular perturbation theory. The analogous case of the 
bacterial growth fronts ,29,] indicates, however, that such 
a derivation can be quite subtle for nonlinear diffusion 
problems. Indeed it is not entirely clear whether the as- 
sumptions used in the sharp interface limit of Ref. [l^lT^ 
are fully justified. For this reason, we have developed an 
alternative and more rigorous stability analysis which al- 
lows for a systematic study on fronts in vortex dynamics. 



B. The model 

The physical situation that we have in mind refers to 
a semi-infinite 2D thin film in which there is an initial 
uniform distribution of vortices due to an external field 
H applied along the z direction. By reversing and in- 
creasing the field, a front of vortices of opposite sign 
penetrates from the edge of the film. We will refer to 
the original vortices as antivortices with density n~ , and 
to the ones penetrating in after the field reversal as vor- 
tices with density n"*". In the region of coexistence of 
vortices and antivortices, annihilation takes place. Vor- 
tices are driven into the interior of the superconducting 
sample by a macroscopic supercurrent J along the y di- 
rection due to the gradient in the density of the internal 
magnetic field. Flux lines then tend to move along the 
direction x transverse to the current under the infiuence 
of the Lorentz force on each vortex (see e.g. 0,0) 

F± = ±-Jx,/.oe„ (1) 

c 

where (/)o is the quantum of magnetic flux associated with 
each Abrikosov vortex. We consider the regime of pure 
flux flow in which pinning can be neglected, while the 
viscous damping then gives rise to a finite vortex mobil- 
ity. We follow a coarse-grained hydrodynamic approach 
in which the fields vary on a scale much larger than the 
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distance between vortices. Since the magnetic flux pen- 
etrates in the form of quantized vortices, the total mag- 
netic field in the interior of the thin film can be expressed 
in a coarse graining procedure through the difference in 
the density of vortices and antivortices, 



B = (r^^ 



Po Gz 



(2) 



The dynamical equations for the fields of vortices and 
antivortices are simply the continuity equations 



dn+ 



= -V • (n+v+) - 



(3) 



where the second term on the right represents the annihi- 
lation between vortices of opposite sign. Note that since 
vortices annihilate in pairs, the total magnetic field is 
conserved in the annihilation process. The annihilation 
terms depends on the recombination coefficient t; a sim- 
ple kinetic gas theory type estimate shows that is of 
order of vS^q , since the cross section of a vortex is of order 
^0, the coherence length IS]. The velocity v can be de- 
termined with the phenomenological formula for the flux 
flow regime: 



= ± J X 



(4) 



where the Hall term has been neglected with good ap- 
proximation for a case of a dirty superconductor . The 
drag coefficient r] is given by the Bardeen-Stephen model 
[sij l and generally depends on the temperature of the 
sample. In this paper we neglect this important coupling 
to the temperature, but we will allow the mobility (the 
inverse of the drag) to be anisotropic. In passing, we also 
note that the above linear relation between the current J 
and the flow velocity is often generalized to a nonlin- 
ear dependence |18| . For simplicity, we do not consider 
this case here, but our method can be extended to such 
situations. 

For a type II superconducting material with a 
Ginzburg-Landau parameter k 3> l/\/2, the magneti- 
zation of the sample can be neglected, so that B « H. 
Then, by using the Maxwell equation (in which the term 
related to the displacement currents has been neglected 
with good approximation) 



J = — V X B, 

47r 



(5) 



together with |5J and and substituting into we 
get: 



dn+ 



(6) 



dt 

— = i?V-(n-V(n--77+)) — . (7) 

where the coefficient D is given by I? = (^{{^ixri). This 
is the system of non-linear differential equations which 



governs the dynamics of the vortex-antivortex front. The 
situation that we will study in our analysis is the follow- 
ing. We consider a front of vortices which propagates 
into the superconducting thin film from the left edge at 
X = —Lx in the positive x direction. At a; = —L^^ we 
impose the boundary condition that the density of vor- 
tices is ramped up linearly in time, n+(— La;,i) = Rt. 
This corresponds to the field going up linearly, just as in 
the Bean critical state j23| . We impose also that far right 
at a; ^ oo, n"*" vanishes while n~ approaches a constant 
value riao ■ Through a rescaling of time and length vari- 
ables, the coefficients of the equations (O and Q can 
be set to unity. In particular, it is convenient to rescale 
the time and length variables according to the following 
transformation: 



a; — > — 



^0 




(8) 



We will henceforth analyze the equations © and O with 
Z) = 1 and T = 1. 

As we already mentioned in the introduction, and as we 
shall see in detail below, the above continuum equations 
have a mathematical singularity at the point where n+ 
vanishes. Of course, in reality there cannot be such a 
true singularity and our continuum coarse-grained model 
breaks down at scales of the order of the London penetra- 
tion depth. In particular, the derivative of the magnetic 
field and thus the current J are not discontinuous with 
respect to the space variable, but they decrease exponen- 
tially in a distance approximately equal to the penetra- 
tion depth. Effects like thermal diffusion, the finite core 
size, and the nonlocal relations which are neglected in 
the London approximation all play a role there, and the 
Ginzburg-Landau equation would provide a more appro- 
priate starting point. Clearly, if the dynamical behavior 
of our continuum model would be very sensitively depen- 
dent on the nature of the singularity, then this would be 
a sign that the physics at this cutoff scale would really 
strongly affect the dynamically relevant long-wavelength 
dynamics. In practice, however, this is not the case. First 
of all, our method to do the linear stability analysis is pre- 
cisely aimed at making sure that the singularities at the 
level of the continuum equations does not mix with the 
behavior or pertubations of the front region. Secondly, 
as we shall see there are no instabilities on scales of the 
order of the microscopic cutoff provided by the London 
penetration depth. 



C. The Method 

In our analysis, we first study a planar front which 
propagates with a steady velocity v along the x direc- 
tion. By considering the propagation of the front in the 
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comoving frame, we get a system of ordinary differential 
equations (ODEs) for the vortex and antivortex density 
fields. The derivation of the uniformly translating so- 
lution is discussed in Section III Al As wc will see, the 
profile that corresponds to the planar front for the den- 
sity of vortices is singular. In particular, in the regime on 
which we will focus, the derivative of the vortex density is 
discontinuous at the point where the field vanishes, while 
in the low-velocity regime there are higher order singu- 
larities. As a consequence of this nonanalytic behavior, 
the numerical integration of the equations has to be done 
with care near the singular point. 

In Section Uni we perform a linear stability analysis of 
the planar solution. A proper ansatz consists here of two 
contributions: a perturbation in the line of the singular 
front and a perturbation of the density field. As we will 
see, the presence of an in-plane anisotropy means that 
the (anti)vortex flow velocity is no longer in the same di- 
rection as the driving force acting on the (anti)vortices. 
Hence, contrary to the isotropic case, we have to con- 
sider a component of the velocity perpendicular to the 
driving force. The viscosity is thus represented by a non- 
diagonal tensor and depends on the angle between the 
direction of propagation of the front and the fast growth 
direction given by the anisotropy. By applying a linear 
stability analysis we get a system of equations for the 
fields representing the perturbation. Through a shooting 
method, and by matching the proper boundary condi- 
tions, we are then able to determine a unique dispersion 
relation for the growth rate of the perturbation. In Sec- 
tion ^3 we treat the case of a stationary front, with a 
velocity v = 0. Contrary to the case of a moving front, 
no singularity in the profiles of the fields is present and 
the analysis can be carried out in the standard way. 



II. THE PLANAR FRONT 

A. The equations and boundary conditions 

In this section, we analyze the planar uniformly trans- 
lating front solutions n+ — Uq {x — vt) , — [x — vt) 
which are the starting point for the linear stability anal- 
ysis in the next section. We refer to the system in a 
comoving frame in which the new coordinate is traveling 
with the velocity v of the front, £^ = x — vt. The tem- 
poral derivative then transforms into dt\x = f^tlc ~ vd(^. 
Since the front is uniformly translating with velocity u, 
the explicit time derivative vanishes. In the comoving 
frame system, we consider ^ to vary in the spatial inter- 
val [— i,+oo]. The equations © and iQ become: 



^ TC' -«o-)-«o%, (9) 



— vortices 

— antivortices 




FIG. 1: Profile of the planar front for the density of the vor- 
tices (n"*") and antivortices for the case v = 1. 



This is a system of two ODEs of second order. Motivated 
by the physical problem we wish to analyze, the relevant 
uniformly translating front solutions obey the following 
boundary conditions at infinity: 



lim Ur. 

« — + 00 ^ 

lim nl 



= 0, 



lim 



lim 



dnX 



0, 



0. 



(11) 



It is important to note that the constant Uoo can actually 
be set to unity: by rescaling the density fields as well as 
space and time, any problem with arbitrary Uao can be 
transformed into a rescaled problem with rioa — 1- The 
stability of fronts therefore does not depend on rioo, and 
in presenting numerical results we always use the freedom 
to set Uoo = 1. 

On the left, the density of vortices n+ increases linearly 
with time with sweeping rate R. After a transient time, 
because of the annihilation process, the field and its 
derivative vanish. The dynamical equation (|5J) for the n+ 
field then yields 



drip 



(12) 



i.e., we recover the well known critical state result |23| 
that in the absence of antivortices the penetrating n"*" 
field varies linearly with slope —v. Requiring that this 
matches the boundary condition n'^{—Lx,t) = Rt for 
large times at ^ = —L then immediately yields that 
R = v^. It can be easily derived that the density of 
antivortices decays with a Gaussian behavior on the left. 
By using indeed the relation (|12|l for large distances and 
substituting it in (^UJ, we get: 



d^ " d^ 



Ae 



(13) 
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Since the analysis of the planar front profiles and of their 
stability is naturally done in the comoving ^ frame, we 
will in practice use a semi-infinite system in the ^ frame, 
and impose as boundary conditions at ^ = — L 



lim n, 







lim Uft 



0, 

const 3> 1. 



hm '^"o 



hm ^ 



0, 



(14) 



Of course, in any calculation we have to make sure that L 
is taken large enough that the profiles have converged 
to their asymptotic shapes. 



B. Singular behavior of the fronts 

Effectively, Eqs. H6I7() and H9I10|I have the form of dif- 
fusion equations whose diffusion coefficient vanishes lin- 
early in the densities and n~ . As already mentioned 
before, it is well known, from e.g. the porous medium 
equation 0) HE Hll , that such behavior induces singular 
behavior at the point where a density field vanishes (see 
e.g. Ref. j32|). Because we are looking at fronts moving 
into the region where n+ = 0, in our case the singularity 
is at the point where the density vanishes. Let us 
choose this point as the origin ^ = 0. Then the relevant 
front solutions have n^{£,) = for all ^ > 0: see Fig. Q] 

Because tIq (0) ^ 0, the prefactor of the highest deriva- 
tive in the n~ equation does not vanish at ^ = 0, and 
hence one might naively think that n~ is nonsingular at 
this point. However, because of the couphng through the 
diffusion terms, this is not so. By integrating Eq. H1U|) 
over an interval centered around = and using that the 
field values Hq and rtp" arc continuous, one immediately 
obtains that 



lim (^_^ 



= 0. 



(15) 



-AC 



Physically, this constraint expresses the continuity of the 
derivative of the coarse-grained magnetic field |(2J). Math- 
ematically, it shows that any singularity in Uq induces 
precisely the same singularity in : to lowest order the 
two singularities cancel. Fig. ^ illustrates this: one can 
clearly discern a jump in the derivative of at the point 
where Uq vanishes with finite slope. 

Before we analyze the nature of the singularity in more 
detail, we note that because of the nonanalytic behavior 
at ^ = 0, it is necessary to analyze the region ^ < 
where ^ separately from the one at ^ > where 
Uq =0. In the latter regions, the equations simplify 
enormously, as the remaining terms in Eq. (|10|) can be 
integrated immediately. Upon imposing the boundary 
conditions Hll|) at infinity, this yields 



driQ 



^ >0. 



(16) 



Let us now analyze of the nature of the singularity 
at ^ = 0. As the effective diffusion coefficient of the 
n+-equation is linear in analogous situations in the 
porous medium equation suggest that the field n"*" van- 
ishes linearly. This motivates us to write for — 1 <C ^ < 



Ate 



.(c), 



(17) 



where fian(C) is the analytic function which obeys 
Eq. H16|) for all ^. Clearly, the continuity condition H15|) 
immediately implies 



At=A^. 



(18) 



If we now substitute the expansion (|17|l with H18|) into 
Eq. lO for tiq we get by comparing terms of the same 
order: 



0{l): A+iv-n-^) = 0, 

O{0 ■■ HA+ ^A^)- 2n-„" - n-„ = 0. 



(19) 



TT-an = '^"-an/'^CIc=0 5 etc. Likewise, if we substitute 
the expansion into Eq. for Uq , we get 



Here 



Oil): 2vA^-2n~^{A+^A,) 



0, 



(20) 



since the term of order unity involving n^^ cancels in view 
of (|16|l . Higher order terms in the expansion determine 
the coefficients A^ and A^, and other terms like Af 
separately, but are not needed here. Together with H16|l . 
the above equations (|19I2UI) immediately yield 



At 



"an(0)= 1/2, 

AV = -17 4- 



1 

16^' 



(21) 



where for convenience we have now put — 1. 

There are two curious features to note about the above 
result. First of all, always vanishes at the point where 
Uq is half of the asymptotic value Uoo at infinity. Sec- 
ondly, note that At is negative for v > 1/4 and pos- 
itive for V < 1/4. Since the vortex density has to 
be positive, we see that these uniformly translating front 
solutions can only be physically relevant for w > 1/4! 

Since the front velocity in this problem is not dynami- 
cally selected but imposed by the ramping rate R = at 
the boundary, we do expect physically realistic solutions 
with V < 1/4 to exist. In fact, it does turn out that in 
this regime the nature of the singularity changes: instead 
of vanishing linearly, nt vanishes with a f-dependent ex- 
ponent. Indeed, if we write for — 1 ^ ^ < [3J| 

- \^nAt+At^+---h (22) 

«o(0 - m^iA- +A-^ + ---)+n-^{0, (23) 
and substitute this into the equations, then, in analogy 
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with the result above, we find 

n-(0)= 1/2, 

At = A^, 
1 

a= — J - 1 > 1, 



(24) 



{v < 1/4), 



while again for ^ > vanishes. A singular behavior 
with exponent depending on the front velocity v is actu- 
ally quite surprising for such an equation [s^l- However, 
one should keep in mind that this behavior is intimately 
connected with the initial condition for the vortices. 
If one starts with a case where n~ does not approach a 
constant asymptotic limit on the far right, but instead 
increases indefinitely, one will obtain solutions where n+ 
vanishes linearly. For this reason, and in order not to 
overburden the analysis with mathematical technicali- 
ties, from here on we will concentrate the analysis on 
the regime v > 1/4. 

Since our study will limit the stability analysis to fronts 
with velocity v > 1/4 in our dimensionless variables, let 
us check how the scale that we consider relates with the 
realistic values of flux flow velocities. By considering re- 
lations IHl) the velocities are measured in units of: 




47ra2 ' 



(25) 



where we have expressed the viscosity 77 in terms of the 
upper critical field Hc2 and the normal state resistivity 
Pn by using j3l|. Furthermore, a is the distance between 
vortices for ^ 00; thus, since Hao = 4>o'^oo ~ 4>vi/o^, 
it follows that Hoq/Hc2 ~ ^o/^^- For the constant r we 
have used the estimate r^^ ~ v^q discussed after Eq. 
We can then rewrite ()25|l as 



47ra4 



Pn 



(26) 



By considering typical values in Gaussian units 
Pn ~ 10~^^s for the resistivity of the material, a co- 
herence length ^0 ~ 2 10~^ cm for high- Tc compounds 
and a magnetic field Hoo ~ 20G (which corresponds to 
a length a ~ 10"'* cm), our velocity is then measured in 
units of u « 1 cm/ s. This velocity scale is much less than 
values found typically in the flux flow regime, since in the 
presence of instabilities fronts of vortices can propagate 
with much higher velocities of order 10"* — lO^cm/s 0. 
Thus, the regime v> 1/4 is indeed the physical relevant 
one. 



structure due to the fact that the annihilation term does 
not effect the difference n+ — n~ . In order to integrate 
the set of equations H9llUf) . it is convenient to consider 
the following transformations in the variables related to 
the sum and difference of the density fields: 



D 
S 



n 



n 



In these variables, the equations become 



dSo 



d ^ dDo 

d~r°^' 

d ^ dDo 



02 



D2 
^0 



(27) 

(28) 
(29) 



By numerically integrating (|28I29() and looking for the 
solutions which satisfy the boundary conditions above, 
we obtained the uniformly translating front solutions. As 
Fig. n illustrates for v = 1, the profile is singular at the 
point where the density of the n+ field vanishes linearly, 
in agreement with the earlier analysis. 

Because of this singularity, the numerical integration 
of the set (|28I29|I is quite nontrivial. In particular, be- 
cause of the discontinuity in the derivative of the field, 
the system (I28I29|I effectively needs to be solved only in 
the interval [— L, 0[, as the matching to the behavior for 
^ > has already been translated into the boundary con- 
ditions H21|l . The first equation can be straightforwardly 
integrated and by combining it with the second, the set 
reduces to 



dDo 
d^ 

dSo 



-vjDo + rioo) 
So 

Soiv{2Do + Ua 



,))^ + iSl-S'oD'o)/2 



vS^ + iv{D'o + UooDo)) 



(30) 

One can easily verify that in this formulation, the expres- 
sion on the right hand side is indefinite at the singular 
point ^ = 0, as both the terms in the numerator and 
denominator vanish. In order to evaluate the expression, 
it is then necessary to perform an expansion of the nu- 
merator and denominator around the critical point val- 
ues S = —D — ncc/2. From such an analysis one can 
then recover the relations (|21|l which we previously ob- 
tained from a straightforward expansion of the original 
equations. Numerically, we integrate the equations by 
starting slightly away from the singular point with the 
help of the results from the analytic expansion. 



III. FRONT PROPAGATION IN THE 
PRESENCE OF ANISOTROPY 



Sum and difference variables 



A. Dynamical equations 



At first glance, the equations look like two coupled sec- 
ond order equations. However, there is more underlying 



As mentioned before, we are interested in the effect 
that an anisotropy in the vortex mobility could have on 
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the stability of the front. In particular, the motivation for 
such an investigation is the experimental evidence that 
an instability for a flux-antiflux front was found in mate- 
rials with an in-plane ab anisotropy, such as for example 

YBa2Cu30r-s El- 

In a material characterised by an in-planc anisotropy, 
the effective viscous drag coefficient depends on the di- 
rection of propagation of the front. More precisely, the 
mobility defined in Q then becomes a non-diagonal ten- 
sor. This leads to a non-zero component of the velocity 
V perpendicular to the driving Lorentz force. We want 
to investigate whether the non-coUinearity between the 
velocity and the force is responsible for an instability of 
the flux-antiflux interface. In the presence of anisotropy, 
the phenomenological formula Q then has to be replaced 
by: 

v = r'F = TR-'(^l l^RF, (31) 

where F is a constant, a represents the anisotropy coef- 
ficient and R is the rotation matrix corresponding to an 
angle 6 between the direction of propagation of the front 
X and the principal axes x' of the sample. The coefficient 
a varies in the range [0,1] with the limiting case of infi- 
nite anisotropy corresponding to a — s- 0. For a = 1 the 
isotropic case is recovered. The matrix f/''^ is given in 
particular by 




The dynamical equations for the fields n+ and n in the 
presence of anisotropy generalize to 




where the length and time variables have been rescaled 
and the elements k and p depend on the angle 9 through 
the formulas: 

cos9sin9(l — a) a cos-^ 9 + sin^ 9 

k= 2"^-, 2—. (34) 

cos^ 9 + a sin 9 cos^ 9 + a sin 9 

Starting from an initially planar profile derived in Sec- 
tion lll AI we want to study the linear stability of the front 
of vortices and antivortices by performing an explicit lin- 
ear stability analysis on Eq. I|33(l . 

B. The linear stability analysis 

As we have already mentioned in earlier sections, our 
linear stability analysis differs from the standard one, due 



to the presence of a singularity. The type of perturba- 
tion that we want to consider should not only involve 
the profile in the region where n+ vanishes, but should 
also in particular involve the geometry of the front. In 
other words, as Fig.[5]illustrates, we want to perturb also 
the location of the singular line at which the density n"*" 
vanishes. As discussed in more detail in '29'|, the proper 
way to implement this idea is to introduce a modulated 
variable 

((^,y,t) =^ + ee'«^+'^'+*"* (35) 

and then to write the densities in terms of this "comov- 
ing" modulated variable. Of course, the proper coordi- 
nate is the real variable Re C. However, when we expand 
the functions in Fourier modes and linearize the dynam- 
ical equations in the amplitude e, each Fourier mode can 
be treated separately. Thus, we can focus on the single 
mode with wavenumber q and amplitude e and then take 
the real part at the end of the calculation. The profiles 
of the fields and n~ are now perturbed by writing 

n+{C,y,t) = n+ iO + e{nt +tn+){C)e'^y+'^'+^^' (36) 
n-iC,y,t) = no-(C) +e(nr +*n2")(C))e*«^+"™ (37) 

where and Hq are simply the planar front profiles de- 
termined before. Note that since we write these solutions 
as a function of the modulated variable C, even the first 
term already implies a modulation of the singular line. 




FIG. 2: Perturbed front profile for the vortex and antivortex 
density field. The fronts propagate in the x direction and has 
a sinusoidal modulation in the y direction 
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Indeed, the standard perturbation ansatz would fail for 
our problem because of the singular behavior of the front. 
The usual ansatz of a stability calculation 

"+(^,2/,i) -no+(0+eK + *«^)(Oe'''^+"*+*"* (38) 

only works if the unperturbed profiles are smooth enough 
and not vanishing in a semi-infinite region. If we impose 
on our corrected linear stability analysis the conditions 



bounded and 



bounded 



Tin 



(39) 

then as e ^ the perturbations can be considered small 
everywhere, plus we allow for a modulation of the singu- 
lar Hue [29j . 

We next linearize the equations H33|) around the uni- 
formly translating solution according to (|36I37(I . We 
obtain a set of 4 linearized ODEs for the variables 
Di, D2, Si, S2, which correspond, respectively, to the real 
and imaginary parts of the difference and sum variables 
introduced in H27() . These equations, which arc reported 
in the Appendix, depend also on the unperturbed pro- 
files Dq, Sq, which are known from the derivation in Sec- 
tion Moreover, there is an explicit dependence on 
the parameters q,uj,n. 

In order to analyze the stability of the front of vor- 
tices and antivortices, the dispersion relation uj{q),n{q) 
must be derived. This can be determined with a shoot- 
ing method: for every wavenumber q there is a unique 
value of the growth rate uj and frequency Q which sat- 
isfies the boundary conditions related to the perturbed 
front. If the growth rate is positive, a small perturbation 
will grow in time, thus leading to an instability. 



C. The shooting method 

The singularity of the front makes the numerical inte- 
gration difficult to handle, as in the case of the planar 
front. In view of the relations (|39|l . the boundary condi- 
tions 



n+ = 0, 



= 0, 



(40) 



have to be imposed for C = 0. These yield the boundary 
conditions for the variables Di, Si, D2, S2 



Si = -Di 



S2 = -D2 



(41) 



Moreover, by substituting these boundary conditions and 
the relations l|21() for the unperturbed fields in the lin- 
earized equations for Di, D2, Si, 82, the followin g re la- 
tions can be derived for C vanishing from the left j35l |: 



dDi 



dD2 



dC 



= LU + qkD2{0), 



= n- qkDi{0) - 2qk 



dDn 



dC 



(42) 
(43) 



An explicit expression for the derivative of the sum of 
the real and imaginary part of the perturbations 5*1 , S2 
can also be derived from the equations reported in the 
Appendix. In particular, these have the following generic 
form 



dSiiO UiiC) 



dC 



2?i(C)' 



dS2 
dC 



A/'2(C) 
2?2(C)' 



(44) 



which is similar in structure to Eqs. (|30|l : A/i,I?i,A/'2,I?2 
depend on ( through the set of functions 



n g dDo dSo n Q n Q 



and on the parameters q,w, fl. 

The equations (|44|l are not defined at the singular 
point. By substituting the boundary conditions given 
by (|41I43|) , both the numerators Afi , N2 and the denom- 
inators T) 1,7)2 vanish. Again, as with (|3U|I . we encounter 
the problem of dealing with the singularity at ^ = 0. 
This difficulty can be overcome in the same way as in 
Section III Bl for the derivation of the planar front pro- 
file. In particular, we can not start the integration at the 
singular point, but we have to start the backwards inte- 
gration at some small distance on the left of C = 0. We 
do so by first obtaining the derivatives of the fields Si 
and ^2 analytically through the expansion of the equa- 
tions H44|l around the critical point. In the limit C ^ 0, 
this yields the following self-consistency condition for the 
derivatives. 



dSi 
~dC 



^fi\o- 



dS2 
dC 



V' 



2IO" 



(45) 



where , N2 , "D'l , 2^2 denote the derivatives of the corre- 
sponding functions evaluated at the singular point. Once 
these are solved and used in the numerics, the integration 
can be carried out smoothly. Because of the singularity 
at the point C = O7 the derivative of the perturbed fields 
are not continuous there and a relationship for the dis- 
continuity in the derivatives can be derived as was the 
case for the unperturbed fields. In particular, the ex- 
pression H15|l is generalized for the perturbed field. This 
implies that the derivative of the total magnetic field is 
again continuous even at the singularity. 

From the equations for the perturbed fields given in 
the Appendix, the boundary conditions at ^ = —L can 
be derived. Just like the unperturbed field for the an- 
tivortex density vanishes on the left with a Gaussian be- 
haviour according to (|13|) , also the perturbations nj~ and 
vanish as a Gaussian, i.e. faster then an exponential. 

Moreover, since the density of vortices increases lin- 
early asymptotically, we can retain in the equations only 
terms which are proportional to the density of vortices 

. From this we get the following equation for the den- 



sity of the perturbation Sn'^ 



int for C ^ — 1 



d^5n^ 



2iqk 



ddnT' 



' pq 5n^ 



pq 



2^ 
dC 



(46) 
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FIG. 3: Dispersion relation u}{q) for different values of 
anisotropy coefficient a and a velocity t; = 1.0. 



FIG. 4: Imaginary part of the growth rate 0.{q) for different 
values of the anisotropy coefficient a, with velocity v = 1.0. 



The solutions of this equation which do not diverge are 
of the form 



5n 



+ Ce^<, X = tqk+^{q^p~k^)), (47) 



where C is an arbitrary constant and k and p represent 
the coefficients of anisotropy defined in Thus, the 

perturbations decay on the left of the film with a decay 
length Co, such that 



7- = 



fc2. 



(48) 



Note that the decay length becomes very large for small 
q; this type of behavior is of course found generically in 
diffusion limited growth models. Technically it means 
that we need to be careful to take large enough systems 
to study the small-g behavior. From the numerical inte- 
gration it was verified that Eqs. l|17|) and describe 
correctly the behavior of Sw^ at large distance. 

Furthermore, since vortices are absent in the positive 
region, we have to impose that the density of the pertur- 
bation related to the n"*" field, and its derivative in space, 
have to vanish there. Similarly we get a second ODE 
with constant coefficients by considering that the den- 
sity of antivortices is constant at large positive distances. 
Taking again = 1, we get, for 3> 1: 



- {v + 2iqk) 



dSr 



dC 



{pq^ +LU + in)Sn- = 0. (49) 



In order to satisfy the boundary condition, we must con- 
sider the solution which vanishes exponentially. The so- 
lution of this equation which does not diverge is of the 
form 



Sn 



Cie^"^ ,Re{X) < 0. 



(50) 



We applied the shooting method in a 4-dimensional space 
defined by the free parameters -Di(O), D2{0),oj and fl, by 



integrating backward in the interval [— L, 0] and then in 
[0, -l-oo[, looking for solutions of the type H47I50|I . 
By matching the solutions to the boundary conditions 



lim riT = — 

C^-L i_ dC 

lim n, — 0, 

C^-l-00 



lim nt = 0, 

(51) 

lim tIt — 0, 

C-+00 ^ 



we then obtain a unique dispersion relation for the real 
part of the growth rate Lu(q). 



D. Results 

Fig. O represents the dispersion relation for an angle 
6* = 7r/4 and different coefficients of anisotropy a. The 
front is always stable, even in the presence of very strong 
anisotropy, for very low values of a. As the anisotropic 
coefficient a is lowered from above, for fixed wavenumber 
q, the growth rate uj{q) increases, but it is always nega- 
tive. For small q a quadratic behavior of uj{q) is found: 



cq , 9 < 1, 



(52) 



where the (negative) coefficient c depends on the 
anisotropy of the sample. In Fig. 21 we have plotted the 
frequency f2 as a function of the wavenumber q. One 
observes from (|35f) that ^1/q is the velocity with which 
the perturbation of the front shifts along the direction 
transverse to the propagation direction. The behavior of 
57 (g) is linear for low wavenumber q and is proportional 
to the non-diagonal element of the mobility tensor fc. 



f2(g) (X kq, q 1. 



(53) 



For an anisotropy coefhcient equal to one the isotropic 
case is recovered and then fi{q) vanishes identically for 
all wavenumbers. As we have already mentioned, the 
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equations that we have used are vahd at scales larger than 
the cutoff represented by the London penetration depth. 
Anyway, since our results clearly show a stability in the 
large q behavior, our model provides a good description 
for the dynamics of the front. In Fig. Owe plot the growth 
rate as a function of for different values of the angle 
9. Linear regression then gives a slope corresponding to 
the constant c in H52|) , which is half the second derivative 
of the growth rate u with respect the wavenumber at 
q — 0. The dependence of c as a function of the angle 9 
is shown in the lower plot. As the angle 9 increases, the 
front becomes more and more stable. This behavior can 
be understood directly from the form of the equations. 
By applying the transformation 

6*^1- 61, 0<6i< 7r/4, (54) 



-10 1 1 1 1 1 1 r 




FIG. 6; Velocity dependence of half the second derivative of 
cj(g) with respect to q evaluated at g = 0. As the velocity 
increases the front becomes more and more stable. 



the elements of the mobility tensor transform into 
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FIG. 5: (a) Plot of i^'{q'^) as a function of the angle 9. (b) For 
a coefficient of anisotropy a = 0.8 and a velocity v = 1.0, the 
results from linear regression for the slope evaluated at g = 0, 
c = dLo/d{q^), are plotted as a function of 9. 



1 

P 



P 



(55) 



By considering the quadratic relation of Lo{q) for small q 
and the fact that the equations are invariant under the 
transformations q = pq and (|55|1 . it is easy to derive 



^{q)\e ^p'^(^{q)U/2~0, 0<p<i, 



(56) 



which proves that the dispersion relation becomes more 
negative as 9 increases. When the direction of propa- 
gation is that of the fast growth direction the isotropic 
case is recovered. In Fig. we show the dependence of 
the coefficient c as a function of the velocity of the front. 
The front is stable for velocities for which Uq vanishes 
linearly (w > 1/4). Furthermore the front becomes more 
stable with increasing v. As one can easily understand 
from the form of the unperturbed front, the vortex 
density profile becomes steeper with increasing the 
velocity. The limit of infinitely large v corresponds to 
the case of a front of vortices propagating in the absence 
of antivortices. Thus, the results confirm the stability of 
the front without an opposing flux of antivortices. 



IV. STATIONARY FRONT 

As we mentioned in the introduction, we have also 
analysed the case of a stationary front, with v = 0. In 
this case it is easy to derive the unperturbed profiles for 
the densities of vortices and antivortices, since they are 
continuous and do not present any singularities. This 
case was previously studied in f|l7|) and treated in terms 
of a sharp interface limit. Equations H9I1U|I in this case 
simplify to: 



" d£, 



'0 "0 



0, 



(57) 
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FIG. 7: Density profiles for vortices and antivortices in tlie 
stationary case {v — 0). The profiles are smootfi and are 
not characterised by singularities, as was the case for fronts 
propagating with finite velocity. 



d _ d _ 



0. 



(58) 



The profiles of vortices and antivortices are symmetric 
in this case, and outside the interfacial zone the density 
fields can be easily derived analytically. By neglecting the 
annihilation term, the profiles of vortices and antivortices 
have a dependence on the coordinate ^ of the type: 



(59) 



where ] — denotes the region where vortices and 

antivortices overlap, N is the density at (±^i) and C a 
constant. The density of vortices and antivortices decays 
with a Gaussian tail, as can easily be calculated from 
equations l(57l and 1551) . For Eq. l|Ffjl . by considering that 
assumes a Gaussian-like dependence, and from the 
form of H59|) . we get the following equation: 



dn^ duQ 
~dr~df 



(60) 



This yields in a self-consistent way a Gaussian behaviour 
for 71+ : 



1+ w ^e~«'^«(^'/^~^^i\ 



(61) 



where A is a constant. The density profiles for vortices 
and anti-vortices are represented in Fig. 7. The stability 
of the front was studied by following a similar procedure 
as for the moving front. Because of the regular profiles, 
the ansatz (|35ll that we have applied for the case of a 
finite velocity is not required. Thus the linear stabil- 
ity analysis can be carried out in the standard way and 
the linearised equations for the perturbation can easily 
be integrated. We do not explain here the procedure 
in detail, since it is a simplified version of the one dis- 
cussed in the previous section. As Fig. |S1 shows, an in- 
stability is found below a critical coefficient of anisotropy 



ac ~ 0.02. These results confirm previous approximate 
calculations 0|, but, as we have already underlined, 
this coefficient would correspond to an extremely high 
in-plane anisotropy which is not found in any type of su- 
perconducting material. We conclude that this model of 
a stationary front in the presence of anisotropy is insuf- 
ficient to explain the turbulent behaviour that has been 
found experimentally at the flux-antiflux boundary. 



V. CONCLUSIONS 

From our analysis it follows that the planar front of 
vortices moving with a sufficiently large velocity w in a 
superconducting thin film is stable even in the presence 
of strong in-plane anisotropy. For stationary fronts, on 
the other hand, our stability analysis confirms the ear- 
lier approximate analysis of T^, confirming that such 
fronts show an instability to a modulated state in the 
limit of very strong anisotropy. From an experimental 
point of view, the critical anisotropy of this instability 
is very high when compared with real values that can be 
found for materials with both tetragonal and orthorhom- 
bic structure [23, even when a non- linear current- 
electric field characteristic is considered • From a the- 
oretical point of view, the behavior in the limit of small 
but finite v is still open as we have not investigated the 
range < u < 1/4 where the profiles have a noninteger 
power law singularity. It could be that the instability 
gradually becomes suppressed as v increases from zero, 
or it could be that the limit w — > is singular, and that 
moving fronts are stable for any nonzero v. Only further 
study can answer this question. 

Our calculations differ markedly from previous work in 
that we focus on moving fronts from the start, where our 
results follow from a straightforward application of linear 
stability analysis to our model. Taken together, these re- 




FIG. 8: Dispersion relation ijj{q) in the case of a stationary 
front. An instability is found for a critical anisotropy coeffi- 
cient Qc « 0.02. 
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suits lead to the conclusion that a model which includes a 
realistic in-plane anisotropy, but which neglects the cou- 
pling with the temperature, cannot explain the formation 
of an instability at a vortex-antivortex boundary for suf- 
ficiently large front velocities. At the same time, our 
calculations show that the issue of the stability of vortex 
fronts is surprisingly subtle and rich. For example, we 
note the fact that for any front velocity, the value Uq at 
the singular line is exactly noo/2 for any v. Is this simply 
a mathematical curiosity or is the absence of instabilities 
related to this unexpected feature through the boundary 
conditions at infinity? Is the presence of a gradient in 
the antivortex distribution far ahead of the front perhaps 
necessary to generate a long- wavelength front instability? 
These are all still open issues, so clearly it is difficult to 
make general statements about the (transient) stability 
of such fronts in less idealized situations. 

One possible interpretation of the results is that when 
one has a finite slab into which vortices penetrate from 
one side, and antivortices from the other side, a station- 
ary modulated front (anihilation zone) forms in the mid- 
dle for extremely large anisotropics. However, a moving 
front never has a true MuUins-Sekerka type instability, 
since a protrusion of the front into the region of antivor- 
tices is always damped as a result of the increased anni- 
hilation. 

The fact that the turbulent behavior at the interface 
between vortices of opposite sign was found in a tem- 
perature window |12|. shows that the coupling with the 
local temperature in the sample has to be considered. 
It appears that it is necessary to include both the heat 
transport and dissipation in the model. Applying an ap- 
propriate stability analysis to such extended models is 
clearly an important issue for the future. 



VI. ACKNOWLEDGMENT 

We are grateful to Hans de Haan for performing time- 
dependent simulations of the one-dimensional vortex dy- 
namics equations, which showed that we were initially 
mistaken about the nature of the induced singularities 
in the n~ field. In addition we have profited from dis- 
cussions with Peter Kes, Rinke Wijngaarden and Gianni 
Blatter. C.B. is supported by the Dutch research founda- 
tion FOM, and would also like to thank the "Fondazione 



A. Delia Riccia" and "Fondazione ing. A. Gini" for ad- 
ditional financial support. M.H. acknowledges financial 
support from The Royal Society. 



APPENDIX A: LINEARIZED EQUATIONS FOR 
THE PERTURBED FRONT 



From the linear stability analysis we get the linearized 
equations for the variables D and S: 



dDn 




nD2+pq'So iDi 
dSo \ f dDi \ f dDo 
dC 



dC 

— qk 



dC 



dc^ 

dDo 
dC 

dC )\dC 




pq^SoD2 = 



dDo 



(A2) 



dS2 
dC 



'So- 



dC^ 

dSo 



" dC^ 
dDo 
dC 



qk 

dDo 
dC 



2So 



dC 

/ dPi d^Dp 
\d^^^ 



dC 



uj{Si + ^] - nS2+pq^Do { Di 



dDo 



dC 



(A3) 



dSi , „ fdDo\ fdDi 
~v—— + 2 



+Di 



dC 

d^Do 



dC 



qk 



dQ J \ dC 
dD2 



D, 



d^Di 
dC 



2Dn 



2'-^D2 



dS2 
dC 



dC dC 
+ pq'^DoD2 
dDo\ ( dD2 , , ^ 



-qk 



2Do 



dC J \ dC 
dDi d'^Do 



fD2 

dC 



- SoSi + DoDi 
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d^C 
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